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Abstract 

The Numerov method for linear second-order differential equations is generalized 
to include equations containing a first derivative term. The method presented has 
the same degree of accuracy as the ordinary Numerov sixth-order method. A general 
scheme of the application to the numerical solution of the Hartree-Fock equations is 
considered. 
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1 Introduction 

The linear second-order differential equations of the type 

y"(x) + g(x)y'(x) + f(x)y(x) =0 (1) 

occur in various fields of physics. Here we shall mean the problem of the numerical self- 
consistent solution of the Hartree-Fock (HF) equations of motion deduced from the Skyrme 
energy functional || describing ground-state properties of atomic nuclei. For spherical nuclei 
these equations can be reduced to the form ([!]), the first derivative term arising due to the 
radial dependence of the nucleon effective mass. Usually the task is solved by the Runge- 
Kutta method. But this method is not the best one for the HF self-consistency procedure 
because it requires the interpolation of functions g(x) and f(x) in Eq. (||) between the 
grid points where function y(x) is not calculated. Another well-known method (see, for 
example, Refs. [§, |l|]), which was proposed by B. V. Numerov in 1923, is free of the pointed 
difficulty but the first derivative term in Eq. ([I]) precludes from its immediate application. 
Several modifications of the Numerov method (NM) were developed 0, in order to include 
equations of the type ([!]) and more general nonlinear equations. Here another generalization 
of the NM is presented which is most suitable for the HF calculations and yields the same 
degree of accuracy as the initial Numerov method. 

2 Generalized linear Numerov method 

Let us introduce notations: y = y(xo), y± = y(x ± h) and analogously for g(x) and f(x), 
where xq is the fixed grid point, h is the step length. Developing quantities y±, y'±, y± in 



powers of h, we obtain 

y + + y„ = 2y + h 2 y^ + ^ + ^ + O(h 8 ), (2) 

y + -y- = 2hy' + ^' + ^ + O(h 7 ), (3) 

y' + + y'_ = 2y' + h 2 y^+^ + O(h e ), (4) 

y' + -y'_ = 2hy>< + + ^yf + 0(h 7 ) , (5) 

y'i + y'L = 2y'> + h^ + ^ + 0(h«), (6) 

y'i-y'L = 2hy'>' + ^ + 0^). (7) 
In addition, Eq. ([I]) yields: 

y'o = -goy'o-foya, (8) 

y'[ = -9+y'+ - f+y+ , (9) 

y- = n // / // • (io) 
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Let us consider Eqs. (PD-flTUP as a system of eight linear equations for eight unknown quan- 
tities: y' , u'q, i/q, y^\ y'±, y'±. Solving these equations and substituting the found quantities 
2/o and y$ in Eq. @ we get 



T y = T +y+ + T_y_ + , (11) 



where 



5h 2 

To = 2 a-— 6 /o, (12) 

T± = a ± ^(10 c g + g + + 9-) + ^b ± f ± , (13) 

a = (l + ^g^^-^g-^+jgg (g+ + g-), (14) 

f>± = f 1 ± ySoJ U T |s T J + (| J 9o3 T , (16) 

C = (^S^O-^-J + fl) 2 ^-' (17) 

l6 

^ {6) = ^(z/f + S^^ + O^ 8 ). (18) 

In more detail this result can be obtained by the following way. Making use of 
Eqs. (D-(P3), we get from Eqs. (0) and (§) 

y'S = 9 - y ~ + f - y -- 9+y ' + - f+y+ -^ + o(k<), (19) 



2h 6 
id _ 2(ff ^ + f Q y ) - g + y' + - f+y+ - g_y'_ - f_y. 



h 2 

-^ + 0^). (20) 

Substituting these equalities into Eqs. (0)^© we obtain 

h 2 

V+-V- = 2hy' + —(g^y'_ + f-y- - g+y' + - f+y+) 

-^# + 0(/> 7 ), (21) 
h 

y'+ + y'- = Zy'o + tt ^-y'- + f-y- - 9+y'+ - f+y+) 

-^ 5) + 0(h e ) , (22) 
h 

y'+-y- = -~(4g y + 4f y Q + g + y' + + f + y + + g-y'_ + f-y-) 

-^ + 0^). (23) 



It is useful to rewrite Eq. ( |2"1~D in the form: 

y'o - 



y+-y- h , 
— 2h — ^3+y+ + f+y+ - g-y- - f-y- 



7h* 

360 



y^+0(h°). 



(24) 



Substitution for this formula into Eqs. 
equations for quantities y' ± 



and ( p3| ) leads to the following system of two 



a + y' + + a_y'_ 
M+'P-y- 



u . 
v , 



where 



u 



h 



h 2 



2h 4 



3 (f-y--f + y + )-^yi 5) + 0(h« 



- h (y + -y- 

2 h 

^9o(y- - y+) - 3 (4/oi/o + + /-?/-) 
+^-9o{f-y- - f+y+) - ^^fo2/o 5) - + c(/i 7 ) 



The solution of the system (|25|) , (|2"6f) is as follows: 

/3 T u ± a T f 



2a 



(25) 
(26) 



(27) 
(28) 

(29) 
(30) 



where a = |(aq_/3_ + a_/3+), that coincides with the definition (II). 

Substituting the formulas (|30"D into right-hand side of Eq. ( p4[) we obtain after some 
algebra with taking into account Eqs. (|27| ) - (|29|) 



y 



S+y+ - S-y- - S y 
2ah 



where 





h 3 


So = 




s± = 






= 





5/i 
12* 



A 
12 



6 



lT-^ T J f± 



(31) 

(32) 
(33) 
(34) 



and 



Finally, substitution for the found solutions y'±, y' Q (Eqs. (|30|), (|3lD ) into Eqs. (§) 
) yields the explicit formulas for the quantities y' ' and in terms of yo, y±. After 
substitution for these formulas into Eq. (0) and a series of lengthy but straightforward alge- 
braic transformations we arrive at the result ( p]) -(18). Omitting the term in Eq. (|TT1) 
we obtain the recurrence three-point formula of the generalized Numerov method for linear 
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second-order differential equations or, for brevity, of the generalized linear NM (GLNM). 
Clearly this method reduces to the ordinary NM if g(x) = in Eq. ([!]). The local truncation 
error of the GLNM, which is contained in the term R^ e \ is one of the same order h e as the 
error of the ordinary NM. 

The formula ( |3lD enables one to calculate the first derivative if the function y(x) is 
known at the grid points. The local truncation error of order h A is determined by the term 
The more precise formula follows immediately from Eq. 



, _ (3 T hg T ) y'^Th (Ag y' + 4f Q y Q + f + y + + T 3R {5) 

V± ~ 3±hg^ ' {6b) 

R® = j^y!P + 0(h 7 ). (36) 

The presence of the derivatives in the right-hand side is a shortcoming of this formula, 
nevertheless Eq. fl35|) is practical for the evaluation of y'{x) at the endpoints of the grid. 



where 



3 Application to the Hartree-Fock calculations 

Consider a general scheme within which the method proposed can be applied to the numerical 
solution of the HF equations. The HF approximation is a basis of numerous microscopic 
physical theories describing the quantum many-body systems. In general formulation, the 
HF method leads to a system of nonlinear integrodifferential equations. We shall consider a 
special case of the HF equations of motion deduced from the energy functional constructed 
on the base of the zero-range Skyrme forces || which are widely used for the description 
of atomic nuclei properties (see, for example, Ref. ||). The variational principle applied to 
the Skyrme energy functional leads in the case of spherical nuclei to the following system of 
equations (in proper units): 

to! ( x ) 

4W - ^Hz'x(x) + 2m 9 (x)[e A - V x {x)]z x {x) = 0, (37) 

Tflq\X) 

where q denotes a sort of nucleon (proton or neutron), the index A stands for the set of 
orbital quantum numbers (including q), x denotes the radial coordinate, z x (x) is the radial 
wave function, e x is the eigenvalue playing the role of single-particle energy, V\(x) is the 
state-dependent mean-field potential, m q (x) is the nucleon effective mass. 

Comparing Eqs. (0) and (j3~7|), we see that they have the same form and would be 
identical if we put 

vn! ^x\ 

y{x) = z x {x), g(x) = q -— , f(x) = 2m q (x)[e x - V x (x)} . (38) 

m q (x) 

The essential difference consists in the following: actually formula (|37]) stands for the system 
of N coupled nonlinear integrodifferential equations because the quantities m q (x) and V x (x) 
are functionals of the densities which depend in turn on the set {z x , z' x } of all the wave 
functions and their derivatives with A 6 {Ai, . . . , Ajv} (see @). In practice, the system 
of equations ( |3"7| ) is solved by making use of some iteration procedure. The convergence is 
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achieved by averaging of the densities calculated on two successive iterations. The description 
and the analysis of the procedure in more detail are outside the scope of the present paper 
(see, e. g., Ref. ||, and references therein where some relevant methods are discussed). Here 
it is important only that on each fixed HF iteration one has to solve the set of N uncoupled 
linear differential equations which have the same form ( p7|) but with already known functions 
m q {x) and V\{x) determined by the results of previous iterations. 

The numerical integration of Eq. (|37|) in this case can be performed by means of 
the GLNM described above. For the sake of simplicity it is convenient to come back to the 
notations of the preceding section taking into account Eqs. (|38|) . Setting 

Y° ut = y /y + , Y° ut = y_/y , (39) 
Yt = Vo/V-, Yl n = y + /y , (40) 

we obtain the following recurrence relations from Eq. ([LTD (omitting the term R^ 6 '): 

yout = T + /(T -T_Y° ut ) , (41) 
Yj n = T_/(T - T+Y™) . (42) 

Note that these transformations of Eq. ([0]) are similar but not identical to ones of the 
renormalized Numerov method developed in Ref. J|. 

In the proposed scheme Eqs. (|4~lD and fl4"2|) are used for the outward and the inward 
integrations, respectively. The outward integration starts from a point near x — 0. The initial 
value of Y° ut in Eq. (|4l|) is calculated using analytic expansions of the regular solutions z\{x) 
about x = 0. To improve the accuracy of calculations it is practical to decrease the step 
length of the grid near this point. The inward integration starts from an outside endpoint 
of the grid where the irregular solutions of Eq. (|3"T| ) are known analytically. The initial 
value of Y™ in Eq. (f42|) is calculated as the ratio of the Whittaker functions for protons and 
of the spherical Hankel ones for neutrons. The eigenvalue e\ is found from the condition 
yout _ \ jyvn some ma tching point Xq = x m . The reasonable choice for this point is (see, 
e. g., Ref. ||) the approximate position of the last extremum of the wave function z\(x). In 
the end of this procedure, which is performed for each z\(x) separately, the function Z\(x) 
is calculated at the grid points using the ratios Y^ ut , Y™, and the normalization condition: 

dxz 2 x (x) = l. (43) 

Finally, the derivatives z' x (x) are calculated employing Eq. ([31]) at the inside points 
of the grid and Eq. (^) at the endpoints. After this, new approximations to the functions 
m q (x), V\(x), which are used in the next HF iteration, are calculated. 

The algorithm, that was briefly outlined above, has been realized in the computer 
code intended for the Skyrme-Hartree-Fock calculations. The convergence and stability 
of the procedure described were tested in the calculations of ground-state properties of all 
doubly magic atomic nuclei using most of the present Skyrme-force parametrizations. It 
was obtained that the algorithm based on the GLNM reproduce the known reference results 
within their accuracy. 
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